home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 19 / CU Amiga Magazine's Super CD-ROM 19 (1998)(EMAP Images)(GB)[!][issue 1998-02].iso / CUCD / Programming / LEDA / source / src / graph_alg / _mwb_matching.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-16  |  7.0 KB  |  294 lines
  1. /*******************************************************************************
  2. +
  3. +  LEDA  3.1c
  4. +
  5. +
  6. +  _mwb_matching.c
  7. +
  8. +
  9. +  Copyright (c) 1994  by  Max-Planck-Institut fuer Informatik
  10. +  Im Stadtwald, 6600 Saarbruecken, FRG     
  11. +  All rights reserved.
  12. + 
  13. *******************************************************************************/
  14.  
  15.  
  16.  
  17. /*******************************************************************************
  18. *                                                                              *
  19. *  MAX_WEIGHT_BIPARTITE_MATCHING  (maximum weight bipartite matching)          *
  20. *                                           *
  21. *  modified 01/94, M.Paul                               *
  22. *                                                                              *
  23. *******************************************************************************/
  24.  
  25.  
  26. #include <LEDA/graph_alg.h>
  27. #include <LEDA/node_pq.h>
  28.  
  29. static list<edge> compute_MWBM( graph& G,
  30.                              const list<node>& A, const list<node>& B,
  31.                              const edge_array<num_type>& weight );
  32.  
  33.  
  34.  
  35. list<edge> MAX_WEIGHT_BIPARTITE_MATCHING( graph& G,
  36.                         const list<node>& A, const list<node>& B,
  37.                         const edge_array<num_type>& weight )
  38. {
  39.   if( A.size()==0 || B.size()==0 ) 
  40.     return 0;
  41.  
  42.   /* change the graph such that A becomes the node set with
  43.      smaller cardinality (if necessary)
  44.    */
  45.   if( B.size() < A.size() ) {
  46.     G.rev();
  47.     list<edge> mwm = compute_MWBM( G, B, A, weight );
  48.     G.rev();
  49.     return mwm;
  50.   }
  51.   else
  52.     return compute_MWBM( G, A, B, weight );
  53. }
  54.  
  55.  
  56.  
  57. /* return the list of possible endpoints of augmenting paths which
  58.    only need a minimal change of the dual function u
  59.  */
  60. static list<node> dijkstra( const graph& G, 
  61.                 const node_list& free,
  62.                             const node_array<num_type>& u,
  63.                             const edge_array<num_type>& weight,
  64.                                   num_type maxval,
  65.                                   node_array<num_type>& dist,
  66.                                   node_array<edge>& pred,
  67.                             const node_array<int>& side)
  68. { 
  69.   node_pq<num_type> PQ(G);
  70.   num_type a_min = maxval,    // minimal distance to an endpoint of 
  71.          b_min = maxval,    // an augmenting path in A resp. B
  72.            dv, dw, c, uv;
  73.   list<node> alist, blist;    // the lists of minimal endpoints
  74.   node v, w;
  75.   edge e;
  76.  
  77.   forall(v,free) {         // initialize distances and PQ
  78.     dist[v] = 0;
  79.     PQ.insert(v,0);
  80.     if( u[v]<=a_min ) {        // compute initial nodelist for A
  81.       if( u[v]<a_min ) {
  82.         alist.clear(); 
  83.         a_min = u[v]; 
  84.       }
  85.       alist.append(v);
  86.     }
  87.   }
  88.  
  89.   while( !PQ.empty() ) { 
  90.     v = PQ.del_min();        // dist[v] is the final distance of v
  91.  
  92.     dv = dist[v];
  93.     uv = u[v];
  94.  
  95.     /* if v is in A, update the nodelist for A
  96.      */
  97.     if( !side[v] && uv+dv<=a_min ) {    
  98.       if( uv+dv<a_min ) {
  99.         alist.clear(); 
  100.         a_min = uv+dv;
  101.       }
  102.       alist.append(v);
  103.     }
  104.  
  105.     /* stop if the actual distance becomes greater than the 
  106.        distance to an endpoint of an augmenting path already seen
  107.      */
  108.     if( a_min<dv || b_min<dv ) 
  109.       break;
  110.           
  111.     /* if v is in B, update the nodelist for B
  112.      */
  113.     if( side[v] && !G.outdeg(v) ) {
  114.       b_min = dv;
  115.       blist.append(v);
  116.     }
  117.  
  118.     /* perform one step of Dijkstra's algorithm
  119.      */
  120.     forall_adj_edges( e, v ) {         
  121.       w = G.target(e);    
  122.       dw = dist[w];
  123.       c = dv + uv + u[w] - weight[e];
  124.  
  125.       if (c < dv) c = dv;  /* rounding errors: u[v]+u[w]-weight[e]
  126.                                                might become negative */
  127.       if( c < dw ) { 
  128.         if( dw == maxval )
  129.            PQ.insert(w,c);
  130.         else
  131.            PQ.decrease_inf(w,c);
  132.         dist[w] = c;
  133.         pred[w] = e;
  134.       }
  135.     }
  136.   }
  137.  
  138.   /* return the appropriate list of nodes
  139.    */
  140.   if( a_min==b_min )
  141.     alist.conc( blist );
  142.   
  143.   return a_min<=b_min ? alist : blist;
  144. }
  145.  
  146.  
  147.  
  148. /* augment matching in one pass along paths of length 1 and 3 
  149. */
  150. static int mwbm_heuristic( graph& G, 
  151.                    const list<node>& A,
  152.                    const edge_array<num_type>& weight, 
  153.                    node_array<num_type>& u)
  154. { 
  155.   node v, w;
  156.   edge e, e2, eb;
  157.   int n = 0;
  158.   num_type max, max2, we;
  159.   node_array<edge> sec_edge(G,0);
  160.  
  161.   forall( v, A ) {            // for all nodes in A ...
  162.     max2 = 0; max = 0; eb = 0;
  163.     forall_adj_edges( e, v ) {        // compute edges with largest
  164.       we = weight[e]-u[target(e)];    // and second largest slack
  165.       if( we >= max2 ) {
  166.         if( we >= max ) {
  167.       max2 = max;  e2 = eb;
  168.           max = we;  eb = e;
  169.         }
  170.     else {
  171.       max2 = we;  e2 = e;
  172.     }
  173.       }
  174.     }
  175.  
  176.     if( eb ) {
  177.       w = target(eb);
  178.       if( !G.outdeg(w) ) {        // match eb and change u[] such
  179.         sec_edge[v] = e2;        // that e2 also gets slack 0
  180.         u[v] = max2;
  181.         u[w] = max-max2;
  182.         G.rev_edge(eb);
  183.         n++;
  184.       }
  185.       else {                // try to augment matching along
  186.         u[v] = max;            // path of length 3 given by sec_edge[]
  187.     e2 = G.first_adj_edge(w);
  188.     e = sec_edge[target(e2)];
  189.     if( e && !G.outdeg(target(e)) ) {
  190.       G.rev_edge(e); G.rev_edge(e2); G.rev_edge(eb);
  191.       n++;
  192.     }
  193.       }
  194.     }
  195.   }
  196.  
  197.   return n;
  198. }
  199.  
  200.  
  201.  
  202. /* compute a maximum weight matching in G
  203.  */
  204. static list<edge> compute_MWBM( graph& G, 
  205.                 const list<node>& A,
  206.                 const list<node>& B,
  207.                 const edge_array<num_type>& weight )
  208. {
  209.   num_type MAX, MIN;
  210.   Max_Value(MAX); Min_Value(MIN);
  211.   
  212.   num_type d_min;
  213.   node v, v_min;
  214.   edge e;
  215.   node_array<edge>     pred(G,nil);
  216.   node_array<num_type> dist(G,MAX), u(G,0);
  217.   node_array<int>      side(G,1), used(G,0);
  218.   int mark=0;            // nodes v with used[v]<mark are unused
  219.   node_list free;
  220.   list<node> vlist;
  221.   mwbm_heuristic(G,A,weight,u);
  222.  
  223.   forall(v,A) {
  224.      if( G.indeg(v)==0 && G.outdeg(v) ) {
  225.        free.append(v);
  226.        dist[v] = 0;
  227.       }
  228.       side[v]=0;
  229.   }
  230.   
  231.   while( !free.empty() ) {
  232.     mark++;            // mark all nodes as unused
  233.  
  234.     /* compute the list of possible endpoints
  235.      */
  236.     vlist = dijkstra(G,free,u,weight,MAX,dist,pred,side);
  237.  
  238.     forall( v_min, vlist ) {
  239.       v=v_min; 
  240.  
  241.       /* if v_min is not the first node of the list, check if the
  242.          augmenting path to v_min contains an used node
  243.        */
  244.       if( v_min!=vlist.head() ) {
  245.         e=pred[v];
  246.         while( e && used[v]<mark ) { 
  247.           v = source(e);
  248.           e = pred[v];
  249.         }
  250.       }
  251.  
  252.       /* if all nodes are unused, augment the matching along the path
  253.        */
  254.       if( used[v]<mark ) {
  255.         v = v_min; e = pred[v];
  256.         while( e ) { 
  257.           used[v]=mark;
  258.           v = source(e);
  259.           G.rev_edge(e);
  260.           e = pred[v];
  261.         }
  262.         free.del(v);
  263.         used[v]=mark;
  264.       }
  265.     }
  266.  
  267.     /* change the dual function
  268.      */
  269.     v_min = vlist.head();
  270.     d_min = side[v_min] ? dist[v_min] : u[v_min]+dist[v_min];
  271.     
  272.     forall(v,A) { 
  273.       if (d_min > dist[v]) 
  274.         u[v] -= d_min-dist[v];
  275.       dist[v] = MAX;
  276.     }
  277.   
  278.     forall(v,B) { 
  279.       if (d_min > dist[v]) 
  280.         u[v] += d_min-dist[v];
  281.       dist[v] = MAX;
  282.     }
  283.   }
  284.  
  285.   /* compute the result
  286.    */
  287.   list<edge> result;
  288.   forall(v,B) 
  289.     forall_adj_edges(e,v) result.append(e);
  290.   forall(e,result) G.rev_edge(e);
  291.   
  292.   return result;
  293. }
  294.